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A numerical method for NNLO calculations 

Gudrun Heinrich*^* 

'^IPPP, Department of Physics, University of Durham, Durham DHl 3LE, England 

A method to isolate the poles of dimensionally regnlated multi-loop integrals and to calculate the pole coeffi¬ 
cients numerically is extended to be applicable to phase space integrals as well. 


1. Introduction 

In the last decades, the interplay between in¬ 
creasing precision from the experimental side and 
NLO predictions becoming available from the the¬ 
oretical side has led to impressive tests of the 
Standard Model. However, certain processes or 
observables in QCD should be known at NNLO 
in perturbation theory in order to match the ex¬ 
perimental precision, and in the domain of elec- 
troweak interactions, many NLO processes still 
await for being calculated exactly. 

The calculation of radiative corrections for elec- 
troweak processes is mainly complicated by the 
presence of several mass scales, whereas in QCD 
the challenge comes from the presence of infrared 
singularities due to massless partons. If we have 
to deal with a combination of strong and elec- 
troweak interactions, i.e. with singularities and 
multi-scale problems, the situation is even more 
involved, and the complexity of the corresponding 
analytic calculations is enormous. On the other 
hand, the performance of computers is improv¬ 
ing continuously, and in order to calculate cross 
sections one generally has to use numerical inte¬ 
gration at some point in any case. All this points 
to the fact that a numerical evaluation of the in¬ 
gredients needed to calculate cross sections be¬ 
yond leading order, like loop integrals and cer¬ 
tain phase space integrals, will be of increasing 
importance in the future. 

In this article I will first present a method 
which has been developed in [1] to evaluate IR (or 
UV) divergent multi-loop integrals by extracting 
the poles in 1/e and calculating the pole coefh- 
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dents numerically. Then I will outline an exten¬ 
sion of this method which makes it applicable to 
the calculation of phase space integrals as well. 
Examples from e+e" ^ 2 jets will illustrate the 
procedure. 


2. The method 


The method, called ’’sector decomposition” [2], 
is a systematic way to disentangle overlapping di¬ 
vergent regions in parameter space. The algo¬ 
rithm which has been developed to automate this 
procedure consists of four basic building blocks 
which will be outlined in the following. 

Consider a scalar graph G with N propagators 
and L D-dimensional loop momenta, typically a 
master integral. The propagators can have pow¬ 
ers Vj > 1 and not necessarily integer. After 
Feynman parametrisation, the graph can be writ¬ 
ten as 
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A necessary condition for the presence of infrared 
divergences is JF = 0, whereas U = Q can only 
lead to an UV divergence. Upon integration over 
Feynman parameters, both IR and UV poles will 
manifest themselves as poles in l/e“ (a < 2L), 
stemming from endpoint singularities {xi = 0) 
of the parameters. However, the singular regions 
are in general overlapping, such that a local sub¬ 
traction procedure for the poles cannot be imme¬ 
diately applied. Our algorithm iteratively sepa¬ 
rates the overlapping regions in parameter space 
by the following steps: 


Part II Iterated sector decomposition 

As after part I the overlapping regions in gen¬ 
eral are not disentangled yet, i.e. T or U still 
vanish if a certain set of parameters goes to zero, 
the decomposition into sectors is iterated until 
Ti and lAi contain a constant term. This iteration 
produces k new subsectors of a given primary sec¬ 
tor I, and in each such subsector one has integrals 
of the form 
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where now all singularities are factored out in the 
bracket, while J^ik and Uik do not lead to singu¬ 
larities anymore. 


Part I Generation of primary sectors 

The integration domain is split into N ’’pri¬ 
mary sectors” by introducing 0-functions: 

N N 

d^x = ^ xY\_(^{xi — Xj)9{xi) 

1=1 j=i 

N 

^ G = Y^Gi. 

1=1 

Then the (5-distribution in Eq. (1) is eliminated 
in such a way that the remaining integrations are 
from 0 to 1, which can be achieved by the sub¬ 
stitution Xj = xitj for j < I, Xj = xi for j = I 
and Xj = xitj-i for j > 1. Because of general 
homogeneity properties of T and U, xi factorises 
like 



U{x)^Ui{t)xf , T{x)^Ti{t)x];+^ (2) 

and thus, using / ^ 5{l-xi{l + Y,k^i 4)) = 1, 
one obtains in each primary sector I 
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Part III Extraction of the poles 

A, < 0 in (3) leads to poles which have to be 
subtracted. Typical for gauge theories are log¬ 
arithmic singularities, i.e. Aj = — 1, where the 
subtractions are of the form 
1 
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= O.e) 
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For Aj < —1, the Taylor expansion around tj = 0 
has to be carried out to a higher order, but the 
procedure works analogously. 

After having isolated the poles in this way, the 
resulting expression can be expanded in e. This 
leads to a Laurent series for each subsector inte¬ 
gral Gik , where the expansion in e can in principle 
be carried out to arbitrary order m: 


2L 

Gik = '^2, ^''=0 ^ 

j = -m 


Finally, the Ri subsectors and the N primary sec¬ 
tors are summed over to obtain the result for the 
graph G: 


N N Ri 
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Note that the feature that singularities only oc¬ 
cur as some Feynman parameters go to zero is 
preserved by the transformations above. 
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Part IV Calculation of the pole coeffi¬ 
cients 

What remains to be done is the calculation 

N Ri 

of the pole coefficients Cj = Yl Yl, ■ These 

are {N — 1 — J)-dimensional parameter integrals 
(j=max(j,0)). In principle, one can attempt to 
perform these integrations analytically, and doing 
this in an automatized way poses no problem for 
large values of j, i.e. the leading and subleading 
poles. However, for smaller values of j, automatic 
analytic integration with standard algebraic inte¬ 
gration routines is not possible anymore. On the 
other hand, numerical integration of these func¬ 
tions does not pose any problem as long as the 
function T has no singularities^ within the inte¬ 
gration region. Therefore, for integrals with more 
than one scale, the numerical points to be calcu¬ 
lated are chosen to be in the Euclidean region in 
order to insure that IF is a regular function. In or¬ 
der to deal with general physical kinematics, more 
sophisticated numerical integration methods have 
to be used, as for example the ones suggested in 
[3,4], which can deal with thresholds in the multi¬ 
leg one-loop case. 


3. Examples for loop integrals 

The method outlined above has been applied 
to check the results for the planar [5] and non- 
planar [6] on-shell massless double box. A pre¬ 
diction [1] for the master integrals of the mass¬ 
less double box with one external leg off-shell has 
also been made, which has been confirmed by an¬ 
alytical calculations later [7-10]. Massless double 
boxes with two off-shell external legs also have 
been calculated. For example, for the planar dou¬ 
ble box with p 3 and P 4 off-shell, one obtains ana¬ 
lytically for the leading and subleading poles: 
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Numerically e.g. at the point 

= (2/3, 2/3, 2/3,1,1), 

the result is 

(2/3, 2/3, 2/3,1,1) = r2(l + e) 

/ 0.8437 2.0524 10.52 48.62 
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The accuracy is better than 3% and can be im¬ 
proved easily at the expense of more integration 
time. 

The most challenging loop integral tackled with 
this method is the planar 3-loop box. It has been 
calculated recently by Smirnov [11] who achieved 
a fully analytical result for the coefficients of 
1/c'^ )/ = 2 .. .6. The numerical result obtained 
with our method is in agreement with [11]. 

4. General multi-parameter integrals 

As the method of sector decomposition is very 
general and straightforward, it can be applied to 
integrals other than loop integrals, as for exam¬ 
ple phase space integrals, as well. However, some 
properties which are specific to loop integrals, and 
thus were built into the code for loop integrals, 
cannot be used anymore in the case of general 
parameter integrals: i) there is not necessarily a 
<5(1 — constraint as in (1), ii) the uni¬ 

versal scaling properties (2) for T and U are lost, 
and - most importantly - in) the singularities 
in e may not only be at Xj = 0 anymore. Tak¬ 
ing into account i) and ii) is not a big issue, al¬ 
though it lengthens the code substantially. On 
the other hand, the absence of singularities other 
than the ones at Xi = 0 is crucial for the program 
to work. However, endpoint singularities can al¬ 
ways be remapped by simple variable transfor¬ 
mations such that they occur at cci = 0 only. As 
in the case of loop integrals, singularities within 
the integration region have to be avoided, for in¬ 
stance by transforming the integrand correspond¬ 
ingly before feeding it into the code. 


- 140.77) . 


^In contrast to the 1/e poles which have been extracted 
already, these are integrable singularities, like for example 
threshold singularities, but they lead to peaks in multi¬ 
parameter space which may cause numerical problems. 


5. Example from e+e —> 2 jets at NNLO 

In order to calculate a cross section like 
e+e” ^ 2 or 3 jets at NNLO, one needs the 
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following ingredients: The two-loop virtual part, 
single radiation from one-loop graphs, and dou¬ 
ble radiation from tree graphs. A lot of progress 
has been achieved in the past 2 years concerning 
the first two items. For the case of double radia¬ 
tion where two particles are unresolved, the soft 
and collinear limits are known, but a systematic 
procedure to set up a local subtraction scheme 
which allows to isolate and analytically integrate 
the infrared singular regions in phase space has 
not been established yet. Here it is shown that 
our code is able to isolate the poles and integrate 
the pole coefficients numerically. The following 
IR singular limits can be distinguished: 

1 . three particles collinear 

2 . two pairs of particles collinear 

3. two particles collinear, one soft 

4. two particles soft 

5. a soft qq pair. 

As a specific example, we will consider the triple 
collinear limit. In this limit, phase space and ma¬ 
trix element factorise in the following way: 

d^{pi,P2,P3,P4) (i4>^^^(pi23,P4) X 

|Mp ^ |MbP X 

■5123 : 


where P 123 and p 4 are defined in analogy to the 
dipole method [12] for NLO calculations as 
P 123 = Pi+P 2 +P 3 - T^P4,P4 = j^P4,y = 

S 123 /S (s = iPl +P2 +P3 +P4f = {Pl23 +P4f), 
Mb is the ‘Born’ 1^2 matrix element with fi¬ 
nal state momenta P 123 , P 4 and {Pgig^q^) is the 
spin-averaged triple collinear q gig 2 q 3 splitting 
function. The triple collinear phase space factor 
can be cast into the following form (d = 4 — 2e) 
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where Zi = PiP4/pi23P4 (* = 1,2,3) and the in¬ 
tegral over M is a rescaled transverse momentum 


integration. In the above parametrisation, the 
Mandelstam variables are given by: 
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found in [13,14]. The aim here is only to demon¬ 
strate how the method operates. So let us fo¬ 
cus first on the Abelian part of the splitting 
q —> gig 2 q 3 , described by the splitting function 
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We want to integrate this function over the triple 
collinear phase space factor given in (4). Using 
parametrisation (5) for the Mandelstam variables, 
we see that the y-integration factorises immedi¬ 
ately, and that the u-integration also factorises 
except if there is an angular dependence through 
523- In this case, however, one has to insure that 
no (integrable) singularity in the analytic plane 
is crossed before feeding the function into the nu¬ 
merical integration routine. This can be achieved 
easily by doing the ^-integration analytically. Us¬ 
ing 

^ [ d9{su4B)~'^’^[\-\-±2 vcosB]~^ 

'\/^r(^ j 0 

= 0(l-i/2)2Fi(l,l + e,l-e,i22) 

-f — 1 )—2 21^1 (1,1-f e, 1 — e, l/i^^) (6) 

one can split the subsequent integrations into two 
parts at = 1. After remapping to integrals 
from 0 to 1 and using the integral representation 
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for the Hypergeometric function, one arrives at 
the ‘standard’ form required by the numerical in¬ 
tegration routine. The result is given by 

= J 

Cl ^ P, 

(47r)ii r^(l — e) e® 

Pi = 4 ± 2 X 10"® 

Pz = 13.99999 ± 4.8 X 10"® 

Pa = 10.1125 ±0.0015 
Pi = -38.8109 ±0.012 
Po = -116.47 ±0.07 

The same procedure can be followed for the other 
singular limits of the matrix element. 

The most complicated denominator occurs in 
the non-Abelian part of the triple collinear split¬ 
ting function, where an angular dependent de¬ 
nominator appears quadratically, 
iPg°g 2 ~qt '^)!^123 ~ ^'(^)/si 2 ± less complicated. 
However, this does not present a problem for the 
numerical method. Again, all one has to do is to 
insure that no singularities in the analytic plane 
are crossed. Carrying out the angular integration 
leads to 

V^r(2 “ e) Jo 

[l±A^±2Acos6»]"^ = 

6>(1-A^) 2Pi(2,2±e,l-e,A2) 

± 0(A^ — 1)^ 2±i( 2,2 ± e, 1 — e, 1/A^) 

Following the same procedure as outlined for (6) 
again leads to a representation which can be fed 
into the numerical routine. 

6. Conclusions 

A constructive algorithm to isolate infrared 
singularities from multi-loop integrals and phase 
space integrals has been presented. The algo¬ 
rithm produces finite parameter integrals as pole 
coefficients, which can be integrated numerically, 
or analytically in simple cases. For the numerical 
integration, one has to insure that the function 


does not have (integrable) singularities within the 
integration region, which can be achieved easily 
by appropriate variable transformations in the 
case of phase space integrals with one overall 
scale. As applications, a result for a massless 
double box with two legs off-shell and an exam¬ 
ple from the phase space of e+e" ^ 2 jets have 
been given. The method can serve to check vari¬ 
ous kinds of analytical NNLO results numerically 
and also provides a step towards a completely nu¬ 
merical evaluation of radiative corrections. 
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